{
    #include "AtlasStyle.C"
    #include "AtlasUtils.h"
    #include "NAssymetricError.C"

    gROOT->SetStyle("Plain");       //no grey background on plots
    gStyle->SetOptStat(0);
    SetAtlasStyle();

    gROOT->LoadMacro("AtlasUtils.C");

    Double_t d0sig[] = {4,5,6,7};

    Double_t d0sig_norm = 5;
    Double_t qcd_norm   = 0;

    // (0: 1 jet incl., 1: 1jet excl., 2: 2jet excl., 3: 3jet excl., 4: 4jet incl.)
    Int_t    ijet       = 4; 

    Int_t n = sizeof(d0sig)/sizeof(Double_t);

    Double_t qcd[n];
    Double_t qcd_err[n];
    Double_t d0sig_err[n];

    bool btagging_cut = false;
    bool triangle_cut = false;

    TString histogram_basedir = "histograms/20110317/";

    ostringstream stm_title;
    stm_title << "sys_err_rctl_j" << ijet;
    if ( btagging_cut ) stm_title << "_btag";
    if ( triangle_cut ) stm_title << "_tria";
    TString plot_title = TString(stm_title.str());

    for (int i=0; i < n; i++) {

        ostringstream stm;
        stm << histogram_basedir << "qcd_mm_d0sig_" << d0sig[i];
        if ( btagging_cut ) stm << "_btag";
        if ( triangle_cut ) stm << "_tria";
        stm << ".root";
        TString fname = TString( stm.str() );

        TFile * f     = new TFile(fname);
        TH1D  * h_qcd = (TH1D*) f->Get("data10_njets_qcd_estat_j0"); h_qcd->SetDirectory(0);
        TAxis * x     = h_qcd->GetXaxis();

        Double_t qcd_4j     = 0.;
        Double_t qcd_4j_err = 0.;
        for(int j=1; j <= h_qcd->GetNbinsX(); j++) {
            Double_t N_jets    = x->GetBinLowEdge(j);
            Double_t N_qcd     = h_qcd->GetBinContent(j);
            Double_t N_qcd_err = h_qcd->GetBinError(j);

            cout << N_jets << " " << N_qcd << " " << N_qcd_err << endl;

            switch(ijet) {
                case 0:
                    qcd_4j     += N_qcd;
                    qcd_4j_err += TMath::Power(N_qcd_err,2);
                    break;

                case 4:
                    if (N_jets >= 4) {
                        qcd_4j     += N_qcd;
                        qcd_4j_err += TMath::Power(N_qcd_err,2);
                    }
                    break;

                default:
                    if (N_jets == ijet) {
                        qcd_4j     += N_qcd;
                        qcd_4j_err += TMath::Power(N_qcd_err,2);
                    }
                    break;
            }
        }

        qcd_4j_err = TMath::Sqrt(qcd_4j_err);

        qcd[i]       = qcd_4j;
        qcd_err[i]   = qcd_4j_err;
        d0sig_err[i] = 0.5;

        if ( d0sig[i] == d0sig_norm ) qcd_norm = qcd[i];

        f->Close();

        delete f;
    }

    // printing the variations of QCD estimation vs. control region
    Double_t stdev = 0.;
    cout << "========================" << endl;
    for( int i=0; i<n; i++) {

        stdev += TMath::Power( qcd[i] - qcd_norm, 2);

        Double_t sys_err = fabs( qcd[i] - qcd_norm ) / qcd_norm;
        cout << d0sig[i] << " " << qcd[i] << " " << qcd_err[i] << " " << sys_err << endl;
    }
    stdev = TMath::Sqrt(stdev);
    Double_t stdev_ratio = stdev * 100.0 / qcd_norm;
    cout << "========================"                          << endl;
    cout << "std. dev: " << stdev << " (" << stdev_ratio << "%)" << endl;
    cout << "========================"                          << endl;

    // creating the plot
    c = new TCanvas("control",plot_title,200,10,700,500);
    //c->SetFillColor(42);
    //c->SetGrid();
    c->GetFrame()->SetFillColor(21);
    c->GetFrame()->SetBorderSize(12);

    // create the TGraphErrors and draw it
    g1 = new TGraphErrors(n,d0sig,qcd,d0sig_err,qcd_err);
    g1->SetTitle("statistical uncertainty");
    g1->GetXaxis()->SetTitle("d0 significance");
    g1->GetYaxis()->SetTitle("number of events");
    g1->GetXaxis()->SetLabelSize(.03);
    g1->GetYaxis()->SetLabelSize(.03);
    g1->GetXaxis()->SetTitleSize(.03);
    g1->GetYaxis()->SetTitleSize(.03);
    g1->GetXaxis()->SetTitleOffset(1.5);
    g1->GetYaxis()->SetTitleOffset(1.5);
    g1->SetFillColor(42);
    g1->Draw("AE3");

    g2 = new TGraphErrors(n,d0sig,qcd,d0sig_err,qcd_err);
    g2->SetTitle("qcd estimation 4 jets incl.");
    g2->SetMarkerColor(kBlack);
    g2->SetMarkerStyle(20);
    g2->Draw("SAME PZ");

    // drawing horizontal lines showing the deviations based on the norminal value
    TLine *l1 =  new TLine(d0sig[0], qcd_norm+stdev, d0sig[n-1], qcd_norm+stdev);
    TLine *l2 =  new TLine(d0sig[0], qcd_norm-stdev, d0sig[n-1], qcd_norm-stdev);

    l1->SetLineWidth(2); l2->SetLineWidth(2);
    l1->SetLineColor(2); l2->SetLineColor(2);
    l1->SetLineStyle(2); l2->SetLineStyle(2);
    l1->Draw("SAME");    l2->Draw("SAME");

    ostringstream stm;
    stm << "#sigma_{d0sig=" << d0sig_norm << "} = " << stdev;
    TString stdev_str = TString( stm.str() );
    myText(0.65,0.8, kRed, stdev_str);

    c->Update();

    c->SaveAs(plot_title + ".eps");
    c->SaveAs(plot_title + ".png");
    c->SaveAs(plot_title + ".pdf");
}
